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We study the accuracy of the post-Newtonian (PN) approximation and its formal region of validity, 
by investigating its optimal asymptotic expansion for the quasi-circular, adiabatic inspiral of a point 
particle into a Schwarzschild black hole. By comparing the PN expansion of the energy flux to 
numerical calculations in the perturbative Teukolsky formalism, we show that (i) the inclusion of 
higher multipoles is necessary to establish the accuracy of high-order PN terms, and (ii) the region 
of validity of PN theory is largest at relative C(l/c 6 ) (3PN order). The latter result suggests that 
the series diverges beyond 3PN order, at least in the extreme mass-ratio limit, probably due to the 
appearance of logarithmic terms in the energy flux. The study presented here is a first formal attempt 
to determine the region of validity of the PN approximation using asymptotic analysis. Therefore, 
it should serve as a template to perform similar studies on other systems, such as comparable-mass 
quasi-circular inspirals computed by high-accuracy numerical relativistic simulations. 
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I. INTRODUCTION 

Before recent breakthroughs in numerical relativity 
(see e. g. [l[ for a review), the post-Newtonian (PN) ex- 
pansion had long been regarded as the best tool to pre- 
dict the evolution of compact binaries. In spite of these 
great numerical advances, the PN approximation is still 
essential for the early inspiral phase, and especially for 
generic systems with arbitrary spins and eccentricities. 
Therefore, this approximation is still indispensable in the 
construction of templates for gravitational-wave detec- 
tion through Earth- and space-based interferometers. 

An accurate knowledge of the gravitational-wave phase 
is crucial for interferometric detection of compact bina- 
ries. In turn, the phasing accuracy crucially depends on 
the convergence (or divergence) properties of the PN ex- 
pansion 0]. For this reason, the structure of the PN 
approximation has been extensively studied with the fol- 
lowing two goals: to determine the region of validity of 
the series; and to improve its accuracy through resumma- 
tion techniques. We shall not discuss the latter problem 
here, but we refer the reader to @, H i, @, 0, B S E3] • 

Early studies of the accuracy of the PN approximation 
focused on head-on collisions and on the quasi-circular 
inspiral of extreme-mass ratio (EMR) compact binaries. 
Simone et al. [ll|, [12] compared the relative 0(l/c 4 )- 
accurate expansion of the energy flux to numerical per- 
turbative calculations. They found that the series con- 
verges slowly for particles falling radially into black holes, 
remaining accurate for v/c < 0.3, where v is the orbital 
velocity and c is the speed of light. Building on previous 
work [l3|, [2 EH El| , Poisson [I?) compared the energy 
flux for quasi-circular EMR inspirals computed with an 
Oil/c 1 ^-accurate PN expansion of black hole perturba- 
tion equations to numerical results. Poisson found that 
the PN series performs poorly for v/c > 0.2, and that 



higher-order terms do not necessarily increase the de- 
tection performance of the series, as measured by the 
fitting-factor with numerical waveforms. Various authors 
argued that the PN series should converge much faster 
for comparable-mass binaries (see e. g. fl2l [IH [l9|). 

Previous studies of the accuracy of the PN expansion 
have been necessarily limited by the lack of accurate (ide- 
ally exact) numerical solutions of the non-linear Einstein 
equations, especially for the comparable-mass case. As 
increasingly accurate numerical evolutions of compact bi- 
naries become available, such convergence studies should 
be revisited, taking into account the improved knowl- 
edge of the "true" numerical evolution of the system. 
Until now, however, a systematic asymptotic analysis of 
the accuracy of the PN series was lacking. The term 
"asymptotic analysis" does not here refer to the study 
of the structure of the series at future or past null in- 
finity. Instead, we mean those techniques applicable to 
asymptotic series, which arise as approximate solutions 
to non- linear partial differential equations. 

In this paper, we perform an asymptotic analysis of the 
accuracy of the PN approximation by investigating the 
quasi-circular inspiral of EMR compact binaries. Asymp- 
totic methods assume nothing about the convergence of 
the series, but only that it derives from the approximate 
solution to a consistent system of differential equations 
(see e. g. [2(| for an introduction). In particular, we shall 
here search for the optimal asymptotic expansion of the 
PN series, i. e. for the truncation order beyond which the 
error in the series becomes larger than expected (that is, 
larger than the next term in the series). The main goal 
of this paper is to determine the approximate region of 
validity of the PN approximation for for quasi-circular 
EMR systems as a function of PN order. 

The approximate region of validity is bounded by the 
region where the true error in the PN series is compara- 
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ble to the PN error estimate. By "true error" we here 
mean the difference between the PN estimate and the 
"exact solution" (i. e. the numerical result), while the 
"PN error estimate" is simply the next order term in the 
series. The orbital velocity will serve as the independent 
variable that labels this region, since this is a coordinate- 
invariant quantity, which for EMR systems is related to 
the angular velocity via to = v 3 /M, with M the total 
mass of the system. The orbital velocity beyond which 
these two errors become comparable marks the region 
outside which one cannot neglect higher-order terms in 
the PN expansion. This is simply because, for larger ve- 
locities, the next order terms in the series are as large as 
the true error in the approximation. 

Such an analysis requires we study the temporal evo- 
lution of the approximate solution to the Einstein equa- 
tions. This could be achieved by investigating the PN 
expansion of several different quantities, such as the en- 
ergy flux or the metric perturbation. The specific choice 
of PN-expanded quantity should not strongly affect the 
region of validity estimates, since all such quantities are 
expanded consistently to the same order. Here we choose 
to work with the energy flux, which is also an observable 
and a coordinate-invariant quantity. 

The consistency of this analysis hinges, of course, on 
the error contained in the "exact" numerical solution. 
For the asymptotic analysis to succeed, this numerical 
error must be smaller than the PN error estimate. For 
example, if we investigate the region of validity of some 
PN quantity to relative C(l/c 4 ), the numerical error 
must be smaller than the terms of relative O(l/o 5 ) in 
the series. The energy flux of EMR systems can be accu- 
rately modeled through the Teukolsky perturbative for- 
malism, which, coupled to Green-function methods and 
spectral integrators, guarantees here a numerical accu- 
racy of 0(1O -6 ). We shall see that this numerical accu- 
racy suffices to perform an asymptotic analysis of the PN 
series to relative ©(l/c 11 ) in the most interesting range 
of the particle's orbital velocity. 
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FIG. 1: Edge of the region of validity for different PN orders. 



The edge of the region of validity of the PN series for 



quasi-circular EMR compact binaries as a function of PN 
order N is given by Fig. [1] These results verify and ex- 
tend those of Simone, et al. [ill [l2l]. The error bars in 
this figure symbolize the uncertainty inherent in the def- 
inition of the edge of the region of validity of any asymp- 
totic series, which is defined more accurately in Sec. HT1 
below. For relative 0(1/ c 6 ) and smaller the PN solution 
seems to have a "convergent" character, since the region 
of validity either increases or remains roughly constant 
with PN order. On the other hand, for larger than rela- 
tive 0(1/ c 6 ) there is a "divergent" behavior, reflected in 
the shrinking of the region of validity with increasing PN 
order. As we shall see, this may be associated with the 
appearance of logarithms at high orders in the PN ap- 
roximation. Similar conclusions were reached by Porter 
S E3] j when studying how to increase the accuracy of 
the PN series through Chebyshev resummation. 

A by-product of this analysis is the determination of 
the minimum number of multipoles needed in the numer- 
ical energy flux to perform any type of comparison with 
the PN approximation as a function of velocity. These 
results are presented in Table [II which provides an easy 
reference to determine how many multipoles are needed 
to have an accuracy comparable to the 0(1/ c) PN ap- 
proximation. For example, if one requires an accuracy 
of relative 0(l/c 6 ), then one need only include up to 
1 = 2 between 0.298 < v/c < 0.408, up to I = 3 between 
0.131 < v/c < 0.298 and up to I = 4 for v/c < 0.131. Re- 
markably, for large velocities one can usually simply look 
at the £ = 2 multipolc, while for small velocities one must 
include more and more multipoles. Such a result is a 
consequence of individual multipolar contributions being 
fairly velocity-independent in the large velocity regime. 

The analysis presented in this paper provides a general, 
gauge-independent and systematic method to study the 
region of validity of the PN approximation for any sys- 
tem. We concentrate on EMR binaries here for simplic- 
ity, but the method can be straightforwardly extended 
to comparable-mass, spinning or eccentric binaries. This 
method is more systematic and general than that of Si- 
mone, et al. [12|, and it is similar in spirit to the "PN 
diagnostic" scheme 0, [2l|, [13, [H| • In this scheme, how- 
ever, comparisons with numerical simulations are carried 
out considering time- independent, conserved quantities 
in a quasi-equilibrium framework. We hope that the 
analysis presented here provides a template for numerical 
relativity groups to test the convergence (or divergence) 
of the PN approximation, which is of great interest both 
to the PN and data analysis communities. 

The remaining of this paper presents more details of 
our methods and results, and it is organized as follows. 
Section |TT] summarizes the basics of asymptotic analy- 
sis, formally defining optimal asymptotic expansions and 
providing a pedagogical example of how to determine the 
region of validity of an asymptotic series. Section [TTT] re- 
views the energy flux of EMR compact binaries in quasi- 
circular orbits as obtained in perturbation theory and PN 
theory. Section ITVl calculates the region of validity of the 
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PN approximation through asymptotic techniques. Sec- in the numerical energy flux to test PN theory. Sec- 
tion [V] concentrates on the number of multipoles needed tion I VII concludes and points to future research. 
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TABLE I: Minimum velocity at which we need to include a minimum of £ m i n multipoles to get an accuracy corresponding to 
the Nth PN term. 



In this work we follow the conventions of Misner, 
Thorne and Wheeler (24|: the metric has signature 
(—,+,+,+), and unless otherwise specified we use ge- 
ometrical units where G = c = 1. The notation O(A) 
stands for terms of relative order A, where A is di- 
mensionless. Therefore, PN remainders are denoted as 
0(v/c) q = 0(1/ c 9 ) = 0(v q ), where q is an integer. 
When we say a quantity is of pth PN order we mean 
that it is accurate to relative C(l/c 2p ). The standard 
asymptotic notation will also be used extensively in this 
paper, and it is defined in the next section. 



II. ASYMPTOTIC SERIES 

In this section we shall review the basic prop erties of 
asymptotic series. We refer the reader to [2fll25j for more 
details. We begin by defining remainders and asymptotic 
series, and continue with a description of their convergent 
and divergent properties. We then define the concept of 
an optimal asymptotic expansion, using which we can 
estimate the region of validity of an asymptotic series. 
We conclude this section with a pedagogical example of 
such a calculation applied to Bessel functions. 



A. Basic definitions and notation 

Consider some partial or ordinary, linear or non-linear 
differential equation, whose solution is f(x). Consider 
also a partial sum of the form 



N 



S N 



(x) = ^ a « ( 3 



(1) 



ri=0 



We here choose such a partial sum for simplicity, but 
in general there could be a controlling factor, such as 
an exponential function, that multiplies the partial sum. 



Let us then define the remainder erm(x) as 

e (N)(x) = f(x) - S N (x) . 



(2) 



The limit as N — > oo of the partial sum, Soo(x), is 
said to be asymptotic to the function f(x) as x — > xq, 
i. e. f(x) ~ Soo(a;) as (x — > xq), if and only if 



i(N)(x) -C (X - Xq) 



N 



(x -> Xq) 



(3) 



for all N. The symbol ~ here is used exclusively to 
mean "asymptotic to" and is never intended to mean 
"approximately," which shall be denoted via the symbol 
«. Similarly, the symbol <C has a specific definition: if 
f(x) g(x) asi-> xq, then 



hm 44=0- 
x— >x g(x) 

Then, f(x) ~ s 00 (x) as x — > x if and only if 



x^x Sn[X) 



1 . 



(4) 



(5) 



There are more formal definitions of an asymptotic series, 
but this one will suffice for the analysis of this paper. 
From this definition, it follows that all Taylor series are 
asymptotic series as well. 

Such a definition of an asymptotic series has a clear 
physical meaning: a power series is asymptotic to some 
function if the remainder after N terms is much smaller 
than the last retained term as x — > xq. In this sense, the 
PN expansion of any quantity is an asymptotic series to 
the exact solution of the Einstein equations as v — ► 0. 
This is so because one expects there to exist a velocity 
region where the remainder of the PN expansion after N 
terms is much smaller than the last retained term in the 
series. Such a remainder, of course, can only be computed 
once an exact or numerical solution is known. 
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B. Convergent and divergent series 

An important consequence of the above definition is 
that asymptotic series need not be convergent. A series 
is convergent if and only if 

Jim e (JV ) = , (6) 

N — >oo 

for all x inside a given radius of convergence R around 
Xo, i- e. for \x — xq\ < R. The radius of convergence can 
be computed via the standard Cauchy ratio test: 

R= lim (7) 

n->oo a n+ i 

where a n is the nth term in the series. 

The convergence requirement of Eq. ^ is much 
stronger than the asymptotic requirement of Eq. ([3]) . The 
main difference is that in the asymptotic definition the 
remainder need not go to zero as N — > oo. Thus, asymp- 
totic series can be divergent. In fact, such series often ap- 
proach the exact or numerical solution much faster than 
any convergent series. However, if one insists on adding 
higher-order terms to the series for some fixed value of x 
(or if one pushes the approximation to x S> xq), eventu- 
ally the series will diverge. Thus, a correct use of asymp- 
totic series forces us to truncate them before the answer 
deviates too much from the exact solution. 

The concept of convergence is said to be absolute, be- 
cause it is an intrinsic property of the coefficients a n and 
it requires no knowledge of the exact solution f(x). On 
the other hand, the concept of asymptoticity is said to be 
relative because it requires knowledge both of the coeffi- 
cients and of the exact or numerical solution. This is the 
reason why analyses of the asymptoticity of the PN series 
were not possible before the recent numerical relativity 
breakthroughs . 



C. Optimal asymptotic expansion 

Suppose that we only know a limited number of terms 
in a (possibly divergent) series. We want to determine 
the optimal number to include in the partial sum to get 
an answer as close as possible to the exact solution, i. c. 
the so-called optimal asymptotic expansion. In asymp- 
totic language, we are looking for the partial sum that 
minimizes the remainder in the approximation. 

More precisely, the procedure is as follows: (i) Choose 
some fixed value of \x — Xq\, so that the series becomes 
a sum of N numerical coefficients. In PN theory, these 
coefficients will in general be functions of the symmetric 
mass ratio ri = mi 7712/ (mi + m^) 1 , of the spins of the 
binary members and of the eccentricity parameters; (ii) 
For this fixed \x — Xo\, search over the individual coef- 
ficients in the series. Typically these terms initially de- 
crease in magnitude, but eventually diverge; (iii) Let the 
first minimum in the sequence occur when n — M such 



that M < N; (iv) The partial sum of all terms in the se- 
ries up to (but not including) the Mth term is the optimal 
asymptotic expansion. The A/th term is also an approx- 
imate measure of the error in the optimal asymptotic 
expansion, because it is asymptotic to the Mth remain- 
der as x — > xo. Thus, the optimal asymptotic expansion 
minimizes the remainder, since if we kept on adding more 
terms the remainder would diverge. 

This method works nicely when we have a large num- 
ber of terms in the series. Unfortunately, it does not 
work very well in PN theory, where usually only a few 
terms are known. Moreover, this method depends on the 
chosen value of \x — Xo \. In PN theory, it is precisely this 
value of \x — xq\ that we wish to determine. Nonetheless, 
we can adapt the method to find the region of validity of 
any PN series, as we show below. 



D. Region of validity 

Let us invert the premise and look for the region of 
validity as a function of the number of terms kept in the 
series. This region can be defined as the region inside 
which the remainder e^) an( i the (N + I)th term in the 
series are of the same order. If the series is asymptotic 
to the exact solution, we then expect sjv to be the best 
approximation to f(x) for all x < X with 

0[e {N) (x)] =O(x~x f +1 . (8) 

The convergent or divergent character of the asymp- 
totic series determines how x behaves as a function of the 
number of terms kept in the series. If the series is con- 
vergent, like a Taylor expansion, then we expect the ac- 
curacy of the series to increase as more terms are added. 
This is so provided x is inside the radius of convergence, 
but as already discussed, this is always the case because 
Eq. ([6]) is stronger than Eq. ([3]). Therefore, it follows 
that x tends to increase, and in fact approach the radius 
of convergence, as the number of terms kept in the series 
increases. On the other hand, if the series is divergent, 
then the opposite is true, namely as more terms are kept 
in the series, x decreases. This is so because higher-order 
terms tend to diverge faster than lower-order ones. 

The order symbol encodes here a certain arbitrariness 
rooted in asymptotic analysis. In other words, one must 
decide a priori how different the right and left-hand side 
of Eq. ([8]) should be before equality is declared. This 
ambiguity means that there is not a precise value for 
the region of validity of an asymptotic approximation. 
Instead, this concept is ambiguous up to the order sym- 
bol. In spite of this ambiguity, we can still arrive at 
important qualitative conclusions through such an anal- 
ysis, provided we quantify the ambiguity with error bars, 
as done in later sections. 
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E. A pedagogical example: The modified Bessel 
function 

Before turning to an analysis of the PN series, in this 
section we shall clarify the definitions of previous sections 
with a relatively simple pedagogical example. Consider 
then the following differential equation 



One recognizes Eq. © as the modified Bessel equation 
with index v = 5 (see e. g. [26]). The modified Bessel 
equation does not have a known closed-form solution 
valid everywhere in its domain, and the numerical solu- 
tion of this equation is called a modified Bessel function: 
y{x) — I${x). The asymptotic expansion of the modified 

Bessel function as x —> oo, say I§ , can be found at any 
given order N. For N = 12, for example, we have: 
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FIG. 2: Plot of the modified Bessel function (solid) and 
its asymptotic expansion of order N — 7 (dotted), N = 9 
(dashed) and N — 11 (dot-dot-dashed). 

Figure [2] shows the exact numerical solution to the 
modified Bessel function (solid) and its N = 7 (dotted), 
N = 9 (dashed) and N = 11 (dot-dot-dashed) asymp- 
totic expansions. The figure clearly shows two important 
features of this asymptotic series: (i) The series can be 
extremely accurate quite far from its singular point of 
expansion x — oo. In principle, there is no a priori rea- 
son to believe that any of these approximations will be 
accurate in the domain plotted, since x = 5 is far from 
x = oo. Nonetheless, all expansions plotted are quite 
close to the numerical answer, even up to s ~ 3; (ii) 
The higher the order N of the asymptotic expansion, the 
smaller its region of validity. Observe that the N = 11 
(N = 7) expansion roughly deviates from the numerical 
answer when x » 3 (x f» 2). As already discussed, this 
behavior is usually associated with divergent asymptotic 



series near the edge of their region of validity. 

As a check of this statement, let us now determine 
the region of validity of this asymptotic expansion. As 
explained earlier, this region will depend on N, the order 
of the asymptotic expansion. From Eq. ([2]) we see that 
the remainder is given by erm = I5 — Jg \ Noting that 
the TVth order term in the asymptotic series is simply 
j( N ) _ j^" 1 ^ the region of validity is defined by the 
relation [Eq. (SJ)] 

o(/ 5 -/W) = e)(/f +1 )-jW). (ii) 

Here it suffices to search for the intersection between 
these two curves. 




1 1 1 1 1 1 1 1 1 1 1 

4.5 5 5.5 6 6.5 7 
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FIG. 3: Plot of the absolute value of the remainder of the 
N = 11 asymptotic expansion of the modified Bessel function 
(solid), the N = 12 term (dotted) and the N = 13 term 
(dashed). The kink in the solid line is due to the use of the 
absolute value operator. 
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TABLE II: Approximate edge of the region of validity (x = 
x) for different asymptotic expansions of the modified Bessel 
function. We also present an approximate measure of the 
error in x, as well as the fractional relative error in ig , that 
is S4 N) = 100(J 5 - 4 JV) )/ / 5, evaluated at x. 



Figure [3] plots the absolute value of the remainder of 
the N = 11 expansion (solid), the N — 12 term (dot- 
ted) and the N = 13 term (dashed). For x > 5.25 the 
exact error in the N = 11 expansion relative to the nu- 
merical answer is larger than the terms neglected in the 
approximation, while for x < 5.25 the opposite is true. 
Therefore, the intersection of the solid and dotted curves 
defines approximately the edge of the region of validity 
for the N = 11 partial sum: x w 5.25. 

This analysis can be repeated for any N and the results 
are presented in Table HH which reproduces and extends 
results in Table 3.1 of [20] - The error bars presented in 
the third column are given by the difference between x 
and the intersection of the (N + 2)-term with A^th-order 
remainder. They serve as a reminder that all quantities 
are asymptotic in nature, and thus, can only be inter- 
preted in an approximate sense. 

Table [II] presents several features that are of interest. 
First, observe that the region of validity decreases as the 
order of the expansion increases (the region of validity is 
{x, oo}, so if x increases, the region of validity decreases). 
As explained earlier, this is an important feature of diver- 
gent asymptotic expansions, due to higher-order terms 
diverging sooner than lower-order ones as we approach 
the edge of the region of validity. Second, observe that 
higher-order partial sums are much more accurate than 
low-order ones, when evaluated at their respective edges 
of validity. For example, the N = 11 expansion has an er- 
ror of only 0.06% relative to the exact numerical answer 
when evaluated at x w 5.3. This is a sensible feature 
of optimal asymptotic expansions, i. e. higher-order ap- 
proximations should be more accurate than lower-order 
ones. 

The region of validity of an asymptotic expansion 
should be interpreted with caution. In fact, the asymp- 
totic expansion of the Bessel function could be used out- 
side its region of validity, as defined here. The risk of 
using the Ath-order approximation beyond x — Sx is that 
of making errors larger than those supposedly contained 
in the approximation. For example, the iV = 11 expan- 
sion of the Bessel function produces a fractional error of 
(0.2, 3.5, 36.5)% relative to the numerical answer when 
evaluated at x = (4.5,4.0,3.5), which is always larger 
than the next-order term in the approximation. There- 



fore, beyond the region of validity as defined here, ne- 
glecting the next order term in the approximation leads 
to larger and larger errors relative to the exact numerical 
solution. 



III. EXTREME-MASS RATIO INSPIRALS 

In this section we shall apply the previously-discussed 
asymptotic tools to the inspiral of a small compact ob- 
ject around a supermassive black hole (an EMR system). 
We shall thus consider a black hole of mass mi = M 
and a much smaller object of mass to 2 = < Af. We 
shall further focus on the adiabatic quasicircular inspi- 
ral phase, where the radiation-reaction timescale is much 
larger than the orbital period. 

Such a system is an excellent testbed for the methods 
discussed above. The PN approximation of several quan- 
tities is known to very high order in the EMR limit. Fur- 
thermore, we can numerically compute the energy flux 
with perturbative techniques that are very accurate for 
all velocities. 

Before turning to the study of the region of validity 
of the PN approximation, in the following sections we 
briefly review the derivation and accuracy of the numer- 
ical solution in black hole perturbation theory and the 
analytical structure of the PN expansion of the flux. We 
also present a simple graphical comparison of the PN re- 
sults with the numerical solution. 



A. Numerical calculation of the energy flux 

An exact solution is necessary if we want to use 
asymptotic analysis to determine the region of validity 
of the PN approximation in the extreme mass-ratio limit 
(henceforth denoted PN-EMR approximation). This "ex- 
act" solution can be found numerically through the use of 
black hole perturbation theory. Let us then express the 
metric as a background (in our case, the Schwarzschild 
metric) plus a perturbation: 

g a p = g { a} + h a p, (12) 

with the perturbation assumed to be small, i. e. \h a p\ -C 

1 3^3 1 • The Einstein equations are then linearized in h a p 
and rewritten in terms of the perturbed Weyl Scalar 
Sfy^t, r, 8, (f>). In turn, the perturbed Weyl scalar can 
be decomposed into multipolar components. Working in 
the frequency domain and using spin-weighted spherical 
harmonics -2^m(0, <^), these components are 

*/ w (w,r) = ^-jdndt e Mt „ 2 y/ m (M)[r 4 <5* 4 (i,r,M)] , 

(13) 

The coefficients of the harmonic decomposition obey an 
inhomogeneous differential equation, the Bardeen-Press- 
Teukolsky equation, where the source term is the har- 
monic decomposition of the stress-energy tensor of the 
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orbiting particle. Due to the symmetries of the source 
term, the multipolar components satisfy the symmetry 



relation 



*2 Tn (r,w) = (-l) / *J_ Tn (r,-a;). 



(14) 



The Bardeen-Press-Teukolsky equation can be solved in 
the adiabatic approximation with Green-function meth- 
ods in the frequency-domain 27], 28| . Our numerical code 
has been described in detail in [29l [30L l3ll. l32| . 

The perturbed Weyl scalar can be used to reconstruct 
the fluxes of energy and momentum at infinity. This is so 
because these fluxes are related to h a p at infinity which 
can be computed from convolution integrals of ^i m (uj, r) 
and the (Fourier-decomposed) spherical-harmonic com- 
ponents of the stress-energy tensor. For circular orbits, 
the multipolar components of the energy flux Fi m = E^ m 
(where |m| < €) satisfy the symmetry property 



F, 



i-m 



F, 



£m •) 



(15) 



so we can limit consideration to multipoles with m > 0. 
Therefore, when we refer to "the relative contribution of 
the (I, m) multipole to the total flux" we are really consid- 
ering twice Fi m . The angular momentum flux is related 
to the energy flux through the particle's orbital frequency 
fi, related to the orbital velocity via v — (Mui) 1 / 3 . We 
have computed the energy flux including multipoles up 
to f max = 8 for a thousand equispaced velocities in the 
range v G [10~ 2 , visco], where the orbital velocity at the 
ISCO visco = 6- 1 / 2 ~ 0.408. 

The Teukolsky formalism has inherent errors, but it 
possesses some distinct advantages over the PN expan- 
sion. The principal advantage is that no restriction on 
the velocities is necessary: the numerical solution is valid 
to all orders in v. For this reason the numerical solution 
is not only compatible with the PN-EMR approach, but 
it is also very convenient to study the latter's region of 
validity. Formally speaking, the PN-EMR approxima- 
tion is a bivariate expansion in two independent param- 
eters: v <C 1, such that the gravitational field is weak 
and all bodies move slowly (PN approximation); and 
fi/M <C 1, that enforces the EMR limit. Therefore, there 
are two independent uncontrolled remainders of relative 
0(N, M), which stands for errors of relative 0(1/ c) N and 
0(n/M) M . The perturbative approach is valid to rela- 
tive C(oo, 0), whereas the PN-EMR expansion (discussed 
below) is only available up to relative 0(11,0). 

The asymptotic tools discussed in this paper hinge on 
controlling the numerical error inherent in the "exact" 
solution. These errors cannot always be modeled analyt- 
ically (see [13, [H| for a detailed discussion). When we 
consider the energy flux for a circular, adiabatic inspi- 
ral, the dominant errors can be roughly divided in two 
groups: 

(1) Truncation errors, due to neglecting high-£ mul- 

tipoles in the multipolar decomposition of Ft m ; 

(2) Discretization errors, due to solving numerically 

the inhomogeneous Teukolsky equation. 



Truncation errors tend to increase with orbital velocity, 
and to a certain extent they can be modeled analytically. 
Poisson [l3| has shown that, for a given £, the luminos- 
ity is dominated by modes with even t + m, while the 
power radiated by modes with odd t + m is suppressed 
by roughly a factor v 2 . Since the total power radiated in 
the £th multipole scales like 



F t ex v^ 1 ^ 



(16) 



the error on the flux due to truncating at some I = 



can be approximated by SFg 



„2£n 



Equation ((T 



provides a simple rule of thumb to determine the number 
of multipoles required to achieve a given fractional trun- 
cation accuracy euy. we must include multipolar compo- 



£(£.) I27J- According 



nents up to € max , where v^ 21 ™^'^ 

to this rule of thumb, including up to ^ max = 8 (as we 
do, unless otherwise stated) yields a truncation accuracy 
better than em ~ 10~ 6 up to v ~ 0.316, while at the 
ISCO we have an accuracy better than cm «2x 10 -5 . 

Discretization error is kept small by use of adaptive 
ordinary differential equation integrators to solve the 
Bardeen-Press-Teukolsky equation, and Gauss-Legendre 
spectral methods to compute convolution integrals [33|. 
We investigated this error by experimenting with toler- 
ance parameters in the integrators and increasing the 
number of points in the spectral methods. We found 
that discretization error at low velocities is very sensi- 
tive to other details of the code, such as the accuracy 
in the numerical inversion yielding r(r»), where r is the 
Schwarzschild radius and r* is the tortoise coordinate. 
At small velocities and large radii we obtain agreement 
with Ref. [28| to a six-digit level, so we estimate the dis- 
cretization error to be roughly of 0(1O -6 ) (and possibly 
smaller at larger velocities). This is confirmed by the 
low- velocity, low-amplitude region of Fig.[S]below, where 
oscillations due to discretization error appear only for 
v < 10~ 2 . These oscillations are in fact of C(10 -7 ) rela- 
tive to the dominant Newtonian flux. 

Before proceeding to the PN expansion of the energy 
flux in the EMR limit, we must mention that both the 
perturbative calculation and the PN expansion neglect 
absorption of radiation by the black hole. Poisson and 
Sasaki [l4[ have shown that this absorption is negligi- 
ble with respect to the energy carried off to infinity, be- 
ing suppressed by a factor of v 8 . Although these terms 
might modify the energy flux ever so slightly, both the nu- 
merical and PN flux considered here consistently neglect 
them. Therefore, these terms cannot affect any conclu- 
sions regarding the region of validity of the PN approxi- 
mation. 



B. PN expansion of the energy flux 

The PN approximation is an expansion in small veloc- 
ities d«1, but to model EMR systems we also expand 
in n/M <C 1. In this PN-EMR approximation the energy 
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flux is given by [16|, [17 



F = -Pkewt 



' N 

E 

.fc=0 



(a* + b k hiv) v k 



(17) 



where N labels the iVth-order partial sum or PN approx- 
imant (known up to N = 11 in the PN-EMR approxima- 
tion), a\ = b\ = bi = &3 = ^4 = ^5 = bi = 0, and 
the remaining coefficients can be found in Eq. (3.1) of 
Ref. [HI]. As explained in the Appendix, the logarithms 
do not affect the applicability of asymptotic methods in 
the velocity regions of interest. Here v is the orbital ve- 
locity, and the Newtonian flux is given by 



5 



(18) 




FIG. 4: Top panel: plot of the absolute magnitude of the 
coefficients of the flux as a function of PN order in a log-linear 
plot. Bottom panel: plot of the sign of the coefficients of the 
flux as a function of PN order. Crosses and circles denote 
polynomial (at) and logarithmic (bk) terms, respectively. 

The structure of the PN-EMR series is worth dis- 
cussing further. The series is composed of two distinct 
types of terms: those that depend only on a certain power 
of the velocity, and those that include a logarithmic de- 
pendence. Since be ^ 0, the logarithmic terms appear 
first at 0(l/c) 6 . Furthermore, as shown graphically in 
Fig. 21 the sequence of coefficients at or bk does not 
present a clear alternating pattern. Finally, the absolute 
magnitude of the coefficients in the series is not of order 
unity, as expected of a convergent series. Instead, their 
magnitude grows drastically with increasing order. Such 
a pathological behavior is not present in the comparable- 
mass limit of the PN expansion (see e. g. [18l Il9| | for a 
discussion) . 



C. Comparison of PN-EMR and numerics 

The energy flux computed both numerically and within 
the PN-EMR framework is plotted in Fig. [5] as a func- 
tion of orbital velocity. Here and below, we shall present 



error bars in all figures only when these errors are visible 
(e. g. if the y-axis scale is of O(10 _1 ), the numerical er- 
ror won't be presented, since error bars of 0(1O~ 6 ) would 
not be visible). Figure O is very similar to Fig. 1 of [13], 
except that here the numerical result is computed includ- 
ing terms up to £ max = 8 and we show all terms up to 
5.5PN order in the flux. 




FIG. 5: Plot of the total energy flux F computed numerically 
(solid) and at different PN orders in the PN-EMR approxi- 
mation. Odd orders are plotted in a lighter gray, while even 
orders are plotted in darker gray. Line styles are as follows: 
F {2 \ F (3) is dash-dash-dot; F w , F (5) is dash-dot-dot; F (6) , 
F m is dotted; F (8 \ F (9) is dashed; F (10) , F (11) is solid. To 
avoid cluttering, the top panel shows F^ up to F {6) and the 
bottom panel shows F^ 7 ' up to F^K 

For low velocities (v < 0.2) all PN approximants (ex- 
cept for the 1PN curve) agree well with the numerical 
result, while in the high- velocity region (v > 0.2) the 
lower-order PN approximants deviate from the numerical 
answer. Such a behavior is reminiscent of that first dis- 
cussed by Poisson [T^. Typically, and roughly speaking, 
the partial sums seem to approach the numerical result 
as we increase the PN order. The question we shall an- 
swer in the next section is whether this approach occurs 
at the expected rate for the given order. 

The accuracy of the PN approximants is better seen 
in Fig. [SI where we plot the moduli of the remainders 
\p( N ) — _p|. Light and dark gray curves correspond to 
odd and even PN orders respectively. We plot here only 
the region v < 0.014, since for smaller velocities the re- 
mainders are smaller than 0(1O~ 6 ), and hence contam- 
inated by numerical error. Observe that the dominant 
even-order remainders decrease monotonically with order 
until e(4) and 6(6) cross at v < 0.257 (dark gray square). 
Similarly, the dominant odd-order remainders also de- 
crease monotonically until C( 3 ) and e/g) cross at v < 0.265 
(light gray square). This behavior is "convergent" in na- 
ture, but the region of convergence seems to be larger 
when odd or even terms are considered separately. Such 
a result derives perhaps from the fact that even and odd 
terms are physically different, coming respectively from a 



FIG. 6: Moduli of the remainders \F (N) - F\ as a function 
of velocity. Light gray curves correspond to odd PN orders, 
while dark gray curves correspond to even orders. Line styles 
are the same as in Fig. [5] The light gray (dark gray) squares 
mark the intersection of e( 4 ) and erg) (e(3) and £(5))- The 
light gray (dark gray) diamonds mark the intersection of £( 4 ) 
and e(5) (e( 7 ) and £(8))- Numerical discretization errors at low 
velocities are of order 10 -6 (see Section IIIA), therefore re- 
mainders smaller than about 10 -6 have no physical meaning. 



sum over multipoles moments and integrals of these over 
the entire propagation history of the waves. 

When all orders are considered simultaneously, the PN 
approximation still presents a convergent-like behavior, 
but in a reduced region. This reduction is due to the 
crossing of 6(4) and £( 5 ) at v ~ 0.135 (dark gray diamond) 
and the crossing of £(7) and £( 8 ) at v 0.125 (light gray 
diamond). Note also that the "region of convergence" 
should not be confused with the "region of validity." By 
the former we mean the region where the (N + l)th-order 
remainder is smaller than the iVth-order one. The latter 
is the region where the asymptotic expansion is valid (the 
remainder is of the same order as the next-order term), 
and it will be discussed in more detail by asymptotic 
techniques in section ITVl below. 

In the large velocity region of Fig. [51 the convergent- 
like behavior is replaced by a divergent one. For example, 
£(4) becomes more accurate than £( 5 ) for v > 0.135. Sim- 



ilarly, £( 7 ) seems more accurate than £( 8 ) for v > 0.125 
(see also Table Hill in the next section), and in fact £(7) is 
strikingly accurate for 0.2 < v < 0.35 -much better than 



5 (8) 



and 



; (9)' 



This suggests that the region of validity of 



lower-order approximations may be actually larger than 
that of higher-order ones. In other words, the region of 
validity seems to shrink with increasing order. Asymp- 
totic techniques are ideal to study the edge of the region 
of validity for high- velocities. 

The convergent /divergent transition is hard to see in 
Fig. [6l so for clarity Fig. [7] plots the remainder as a 
function of PN order N for four selected values of v. 
Until now, we have always included multipoles up to 
I = 8, but here we explore the multipolar dependence, 
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FIG. 7: Plot of the absolute value of the remainder of the en- 
ergy flux at different values of v versus PN order in a log-linear 
plot. We present results obtained when summing multipoles 
up to I = 2 (circles), I — 3 (squares), i — 4 (diamonds), £ — 5 
(triangles up), I = 6 mode (triangles down) and I — 8 (stars). 



truncating the sum at I = 2 (circles), I = 3 (squares), 
£ = 4 (diamonds), £ = 5 (triangles up), I = 6 (triangles 
down) and I — 8 (star). Observe that for small velocities 
(i. e. v = 0.01 or v = 0.1) the remainder generally de- 
creases with PN order when using up to £ = 8 harmonics. 
This is precisely the convergent- like behavior alluded to 
earlier. Also observe that an increasing number of multi- 
pole moments must be included as a function of PN order 
to perform any type of meaningful comparison between 
PN theory and numerical results. We shall study how 
this number of multipoles depends on velocity in Sec. [V] 



IV. REGION OF VALIDITY 

Before applying asymptotic techniques to determine 
the region of validity, let us show how a quick but incor- 
rect estimate can be made directly through Fig. [6] Let us 
then forget for the moment that higher-order approxima- 
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tions should be more accurate than lower-order ones, and 
simply draw a horizontal line in Fig. [S] at some maximum 
error threshold, such as 5 max — 0.001. The required ac- 
curacy, of course, depends on the application one has in 
mind. One can then guarantee that inside some region 
v < v the error in the ./Vth-order PN expansion is less 
than <5 max . 



N 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


V 


0.04 


0.11 


0.14 


0.14 


0.19 


0.30 


0.22 


0.25 


0.28 


0.29 


SF iN) [%] 


0.10 


0.10 


0.10 


0.10 


0.11 


0.11 


0.11 


0.11 


0.11 


0.11 



TABLE III: For v < v (second row) the error of using the iVth 
order PN expansion (where N is listed in the first row) is less 
than 0.001, i. e. the relative fractional error SF^ ■ ' ~ 0.1% 
(second row). 



We present v in Table IIII1 together with the rela- 
tive fractional error in the energy flux, = (F — 
F( N 1)/F. This table shows, as expected, that the frac- 
tional relative error in the flux is approximately always 
of the same order, 0.1%. This scheme, however, forces 
the llth-order PN approximant to be as accurate as the 
2nd-order PN approximant, which is clearly inconsistent 
with perturbation theory. 

The restriction that higher-order approximations be 
more accurate than lower-order ones is naturally incorpo- 
rated in the asymptotic estimates of the region of validity, 
which we shall perform next. The condition defining the 
edge of the region of validity is 



Qtp _ p{N)\ _ Qrp(N+l) _ p(N)\ 



(19) 



In other words, we have to determine the v for which 
the remainder ceases to be comparable to the next order 
term. 

The remainder and the next-order term possess two 
distinct and generic features when compared to each 
other. Figure [5] presents the n = 2 (n — 6) order re- 
mainders, as well as the n = 3 (n = 7) order term in the 
top (bottom) panel as a function of velocity. The two dis- 
tinct behaviors alluded to earlier are the following: either 
the remainder and the next-order term are of the same 
order for sufficiently low velocities, until eventually the 
curves separate for larger velocities; or the next-order 
term starts off smaller than the remainder, but even- 
tually the curves cross and separate. When the curves 



cross (in the bottom panel, this crossing occurs roughly 
at v s=! 0.275) one can apply the same techniques of the 
previous section and simply define v as the velocity where 
the curves intersect. However, when the curves do not 
cross (top panel), we must extend the methods we used 
for Bcsscl functions. 

The asymptotic analysis definition of the edge of the 
region of validity is inherently ambiguous, depending on 
the precise meaning of the order symbol. Let us then 

10° i""- 




FIG. 8: Top: Log-linear plot of the absolute value of the 
remainder of the N = 2 PN flux in the extreme-mass ratio 
limit (solid) and the N = 3 term (dotted). Bottom: Same as 
top, but for the N = 6 remainder and the N — 7 term. All 
curves have been factored by the leading-order, Newtonian 
expression for the energy flux. 



replace Eq. (19]) by 



F (N+1) _ p(N) 



< 8, 



(20) 



where 8 is some tolerance. One would expect this toler- 
ance to decrease with PN order, since higher-order ap- 
proximations should be more sensitive to smaller differ- 
ences. Let us, however, forget this for the moment and 
demand a constant tolerance S = 10~ 3 . This procedure 
is somewhat similar to picking a maximum error thresh- 
old in the remainders of the approximation, which we 
already explored in Table ITTT1 Here, however, we do not 
demand a constant maximum error, but instead we arbi- 
trarily choose a constant relative difference between the 
remainder and the next term in the series. 



J 



N 


2 


3 


4 


5 


6 


7 


8 


9 


10 


v(8v) 


0.107(0.017) 


0.138(0.021) 


0.140(0.017) 


0.190(0.020) 


0.266(0.045) 


0.222(0.019) 


0.248(0.019) 


0.281(0.018) 


0.292(0.019) 


8F[%] 


1.51 


0.29 


0.11 


0.61 


1.03 


0.02 


0.28 


0.35 


0.16 



TABLE IV: Approximate values of the edge of the region of validity for different orders of the EMRI-PN energy flux. The first 
row lists v and (in parentheses) 5v; the second row lists SF^* 1 = (F — F^ N ^)/F, evaluated at v. 
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Applying this procedure, one obtains the v and its er- 
ror Sv presented in Table HVl Here Sv symbolizes varia- 
tions of the tolerance in the interval 5 = {10~ 3 /5, 5 x 
10~ 3 }, i. e. we evaluate v\ (v 2 ) with Si = 10~ 3 /5 
(S 2 = 5 x 10 -3 ) and define the error as Sv := \vi — v 2 \j2. 
As one can see from Tabic IIV1 the relative fractional er- 
ror decreases on average with PN order, although not 
monotonically. 




FIG. 9: Tolerance as a function of PN order. 

Higher-order approximations should be sensitive to a 
smaller tolerance, which implies that this quantity can- 



not be set arbitrarily. Instead, 8 should be given by the 
error in the difference between the Ath remainder and 
the (A + l)th-order term. This error is presumably of 
the order of the error in the (A + l)th-order term, and 
it can be estimated by the (A + 2)th-order term. This 
quantity, of course, depends on v, but its order can be 
roughly given by its absolute value evaluated in the mid- 
dle of its range v « 0.2. By this method we estimate 
that the tolerance behaves with PN order as shown in 
Fig. [5] Since this estimate of the tolerance is not exact, 
we shall allow S to vary between 5 times and 1/5 times 
its value. In this way, we can determine the sensitivity 
of our results to the choice of 5 and provide error bars as 
done previously. 

The edge of the region of validity can now be com- 
puted using the tolerance criterion defined above. Ta- 
ble [V] presents v, together with an averaged error bar Sv, 
which represents variations of {(5/5.0,5.0 5}. The first 
line in the first row of this table uses the tolerance pre- 
sented in Fig.[5J Observe that the relative fractional error 
in the flux, the first line in the second row of Table [Vj 
decreases on average as the order of the approximation in- 
creases. However, the region of validity seems to increase 
between 2nd and 6th PN order (except as one goes from 
4 to 5), while it seems to decrease between 6th and 10th 
PN order. We shall analyze this behavior in more detail 
shortly. 



N 


2 


3 


4 


5 


6 


7 


8 


9 


10 


v(Sv) 


0.179(0.029) 


0.200(0.014) 


0.207(0.025) 


0.199(0.021) 


0.288(0.007) 


0.207(0.018) 


0.204(0.016) 


0.178(0.010) 


0.160(0.006) 




0.160(0.026) 


0.200(0.014) 


0.214(0.026) 


0.197(0.021) 


0.346(0.019) 


0.214(0.018) 


0.207(0.016) 


0.158(0.005) 


0.179(0.029) 




0.143(0.023) 


0.200(0.014) 


0.221(0.027) 


0.195(0.020) 


0.390(0.024) 


0.222(0.019) 


0.211(0.016) 


0.160(0.101) 


0.160(0.101) 


SF[%] 


6.7607 


-1.3250 


-0.5562 


0.7888 


-1.6931 


-0.0134 


0.0518 


-0.0043 


0.0001 




4.8505 


-1.3250 


-0.6269 


0.7531 


-5.2159 


-0.0151 


0.0592 


-0.0014 


0.0005 




3.5011 


-1.3250 


-0.7065 


0.7103 


-9.6392 


-0.0168 


0.0686 


-0.0015 


0.0001 



TABLE V: First row: edge of the region of validity v and (in parentheses) estimated error Sv. Second row: 5F (N) [%} with the 
tolerance set by the (N + 2)th-order term in the approximation. Top to bottom, we list values corresponding to three successive 
iterations of our attempt to estimate the region of validity (see text). 



Before discussing the implications of these results, let 
us attempt to determine how reliable they are by exper- 
imenting with the choice of S. Let us then continue to 
define the tolerance through the (A + 2)th-order term, 
but this time let us evaluate it at the value of v found 
in the first line of the first row of Tabic fVl Doing so, we 
obtain the values of v presented in the second line of the 



first row of this table. If we iterate this algorithm once 
more and use the second line in the first row of Table W\ to 
evaluate the (A + 2)th-order term, we obtain yet another 
v, given in the third line of the first row of this table. 

The edge of the region of validity seems to be rather 
insensitive to the choice of S, provided the (A + 2)th- 
order term is not evaluated too far from the mean of the 
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domain. Indeed, Table [V] shows that for most values of 
N the value of v given in the first row remains within 
the error bars of the first line of the first row. This is 
not the case for N = (9, 10), because then the values of 
8F obtained are at the level of the numerical accuracy 
of our simulations: SF^ 9 ' 10 " 1 « 1CP 6 . This is also not 
the case for N = 6, because this is an inflection point 
in the behavior of the edge of the region of validity, and 
thus the first iteration forces v outside of this region. We 
therefore conclude that the first line in the first row of 
Table IVl suffices as a reliable, approximate measure of the 
edge of the region of validity of the PN approximation. 

A note of warning is due at this point: the iterative 
scheme presented here should not be confused with a con- 
vergence scheme. One might be tempted to expect that 
as we continue to iterate, v will tend to some definite 
value that will inequivocably define the region of validity 
exactly. This concept, however, is inherently flawed be- 
cause asymptotic series by definition do not possess an 
exact region of validity. As already explained, asymp- 
totic series and their region of validity are defined via 
order symbols, and thus can only be interpreted as ap- 
proximate concepts. 



locity range, until the curves cross around v — 0.265. 
On the other hand, e (7) w F^ - F^ only for v < 0.2, 
and soon after the curves separate. We see then that the 
absolute magnitude of the remainder itself does not de- 
termine the region of validity. As we stress once more, 
this region is defined by the requirement that the error in 
the approximation is smaller than (or of the same order 
as) the next-order term in the series. 

The behavior of the edge of the region of validity 
with PN order was already presented in the Introduc- 
tion (Fig. [T|). Whether there is truly a negative slope 
for N > 9 is difficult to establish due to the truncation 
error, but our results seem to support this hypothesis. 
Nonetheless, we do observe in the figure that for N < 6 
the region of validity increases (as expected of a conver- 
gent series), while for N > 6 it decreases (as expected of 
a divergent series) . Such a result seems to agree with the 
statement that the logarithmic terms in the PN approx- 
imation somewhat destroy the convergence properties of 
the series H, G3- Nonetheless, as we have seen in this 
paper, a divergent series can sometimes be even more 
useful and accurate than a convergent one, provided it is 
used within its region of validity. 
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FIG. 10: Plot of e (6) (solid black), e (7) (dashed black), the 7th- 
order term (solid light gray) and the 8th-order term (dashed 
light gray). The square shows where £(6) and (F' 7 ' — F' 6 ') 
cross. 

Let us now discuss the behavior of the edge of the re- 
gion of validity with PN order in more detail. As we 
already mentioned, there seems to be an inflection point 
at N — 6: on average, v increases for N < 6, while it de- 
creases for N > 6. Such a behavior then forces the N = 6 
approximation to have the largest region of validity. By 
naively inspecting Fig. [6] one could be tempted to con- 
clude that this result is paradoxical, since €(7) < £(6,8)- 
However, recall that the region of validity is determined 
by the difference between erm and the (N + l)th-order 
term (not e(jv+i))- To stress this difference, in Fig.flOlwc 
show £( 6 ), 6(7) (black) and the 7th- and 8th-order terms 
(light gray). Although the absolute value of em is smaller 
than e( 6 ), observe that C(q) ~ F^ — F^ 6 ) in a large ve- 



FIG. 11: Minimum orbital radius [Eq. (|21[) ] representation of 
the edge of the region of validity for different PN orders. 

The edge of the region of validity can also be pre- 
sented as a function of the Schwarzschild radius of the 
particle's circular orbit. Figure 1111 shows the minimum 
Schwarzschild radius delimiting the region of validity, as 
a function of PN order. This radius can be obtained from 
v through the generalization of Kepler's law to circular 
geodesies: 

i = h- (21) 

Simonc et al. determined that the PN approximation 
for the infall of a particle into a Schwarzschild black hole 
can be trusted down to a harmonic radius /"harm ~ 10M 
[ill ]. Perhaps it is no coincidence that their final result 
is surprisingly close to the edge of the region of validity 
for the 3PN approximant of the flux for quasi-circular 
inspirals. 
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The results presented so far should not be misinter- 
preted as saying, for example, that the 3PN approxima- 
tion is more accurate than higher-order ones. In fact, 
the statements made here say nothing about the abso- 
lute accuracy of the PN approximation. These results 
only suggest relational statements between the Nth and 
the (N + l)th-order approximations. Table fVl should be 
read as saying that, for velocities v < v, the Nth order 
approximation has errors that are of expected relative 
size. For larger velocities, the (iV+l)th- and higher-order 
terms become important, and should not be formally ne- 
glected. However, if one is willing to tolerate errors larger 
than those estimated by the (N + l)th order term, and 
if one is willing to relinquish the desire that higher-order 
approximation be more accurate than lower-order ones, 
then one can surely go beyond v at the cost of losing 
analytic control of the magnitude of the error. 



V. MULTIPOLES AND PN ORDERS 

As a by-product of this analysis, we can ask the fol- 
lowing interesting question: how many multipoles should 
be kept in the numerical solution if we are interested 
in studying the Ath-order PN approximation? In other 
words, we wish to determine whether it is sufficient to 
keep only the quadrupole in certain cases, or if we al- 
ways need higher multipoles to study the accuracy of the 
PN series (see [34] for a related discussion). 
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FIG. 12: Plot of the energy flux decomposed into spherical 
harmonics and summed over m. Each multipolar component 
is normalized to the Newtonian flux. The 1 = 2 harmonic 
is shown with a solid line, the t = 3 with a dotted line, the 
I = 4 with a dashed line, the £ = 5 with a dot-dashed line, 
the I = 6 with a dot-dot-dashed line, and the 1 = 7 with a 
dash-dash-dotted line. 

The answer to this question depends on the behavior 
of each multipolar contribution to the energy flux as a 
function of velocity. In Fig. [12] we plot multipolar contri- 
butions with different €s (where for each t we sum over 
all values of m) normalized to the Newtonian flux. The 
dominant (I = 2) mode is very close to unity for all val- 



ues of v, while all other contributions go to zero as v — > 
and approach a roughly constant, ^-dependent value. We 
do not show in the figure the I = 8 harmonic to avoid 
cluttering, but this contribution would be located at the 
lower right corner of the figure. 

Figure [12] implies that it is ludicrous to compare PN 
with numerical relativity in certain systems if only a 
few multipoles are taken into account. For example, ne- 
glecting the t = 4 multipole leads to errors of 0(1O) -4 
at v = 0.1, which is comparable to terms of 0(l/c) 4 . 
Thus, comparing PN expansions of C(l/c) 5 and smaller 
to numerical solutions neglecting the I — 4 multipole at 
v < 0.1 is risky and may lead to incorrect conclusions. 
The relative importance of higher-i' multipoles decreases 
for equal mass systems, but again recovers great impor- 
tance when eccentricity and/or spins are included. 

The number of multipoles to retain also depends on the 
PN contribution to the energy flux. Figure [TBI plots these 
contributions, normalized to the Newtonian flux. As be- 




FIG. 13: Plot of the relative PN contributions to the total en- 
ergy flux normalized to the Newtonian flux. Top: PN orders 
C(l/c 2 )toO(l/c 6 ). Bottom: PN orders C(l/c 7 ) to C(l/c n ). 

fore, odd (even) orders are plotted in light (dark) gray. 
Observe that the PN contributions of the top panel (lower 
PN orders) rise rather fast and approach a roughly con- 
stant value. The opposite is true for the bottom panel, 
where the curves are monotonically increasing in this ve- 
locity domain. 

An analysis of the accuracy of the Ath-order PN ap- 
proximation requires that the numerical error be at the 
very least smaller than the Ath-order term. In the 
asymptotic analysis of the previous section, we simply 
summed up to I = 8 and always included all harmonics 
in our analysis. However, it is possible that we did not 
need to keep up to I — 8 in the analysis of all PN or- 
ders. Figure [JJ] superimposes the PN contribution and 
the multipolar contribution of the energy flux. We see, 
for example, that the I = 5 harmonic contribution is al- 
ways negligible if one is studying a 1.5PN order accurate 
expression. 

The intersection of the harmonic and PN contribu- 
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FIG. 14: Same as Figs. [T21and fTSl The black solid curves are 
the harmonic contributions to the energy flux, while the non- 
solid curves represent the PN contributions. The intersections 
in this figure (some of which are shown as a square, a diamond 
or a circle) lead to Table U in the Introduction. 

tions provide a minimum requirement of accuracy if any 
type of comparison between numerical simulations and 
PN theory is to be carried out. These intersections form 
the basis of Table U in the Introduction. For example, 
if one wishes to compare numerical results to the 1.5PN 
energy flux, then one can simply use the £ — 2 harmonic 
provided v > 0.101 (shown as a circle in the figure). For 
smaller velocities, however, the £ = 3 harmonic contribu- 
tion needs to be included. Similarly, if one is studying the 
2.5PN expression of the energy flux, then one can simply 
use the £ = 2 harmonic if v > 0.298 (shown as a diamond 
in the figure), but one must include the £ — 3 harmonic 
if 0.057 < v < 0.298 and the i = 4 harmonic if v < 0.057 
(a square in the figure). Finally, note that the num- 
ber of multipoles that need to be included is larger for 
smaller velocities. This is because the multipolar contri- 
bution raises steeply for small velocities and then roughly 
asymptotes to a constant, leading to more intersections 
when v is small. 



VI. CONCLUSIONS AND OUTLOOK 

We have proposed a generic, gauge-independent and 
systematic method to determine the formal region of va- 
lidity of the PN approximation. This method relies heav- 
ily on tools from asymptotic analysis, and in particular, 
on the concept of an optimal asymptotic expansion. The 
main philosophy of the approach is to determine the ve- 
locity or frequency beyond which the next-order term 
in the series must be included. This velocity is found 
by studying the region where the true error in the ap- 
proximation (relative to the exact or numerical answer) 
becomes comparable to the series truncation error (due 
to neglecting higher-order terms in the series) . 

We applied this method to the quasi-circular inspiral 
of a small compact object into a non-spinning black hole. 



The PN approximation is thus coupled to an EMR ex- 
pansion, leading to what we here called the PN-EMR ap- 
proximation. The exact solution is modeled by a numer- 
ical calculation in black hole perturbation theory. This 
scheme linearizes the Einstein equations in the metric 
perturbation but not in the velocity parameter, and so it 
is suitable for the study of EMR systems. Here we are 
mainly interested in a proof-of-principle of the proposed 
method and we used the energy flux as the dependent 
variable, but we could have studied other time-dependent 
quantities. The gravitational wave phase, for example, is 
the most interesting quantity for data analysis purposes. 
However, relating the energy flux to the phasing involves 
a certain degree of arbitrariness (see [35| for a discus- 
sion) and we preferred to avoid these complications in 
our preliminary exploration. 

The output of the method is the edge of the region of 
validity as a function of PN order (Fig. [T|) . We found that 
the 3PN approximation of the energy flux seems to pos- 
sess the largest region of validity. This implies that the 
3PN approximation can be evaluated up to rather high 
velocities (i. e. v w 0.29) and still produce an approxi- 
mate answer with an error of expected size (given by the 
next order term in the series). Moreover, we found that 
for lower than 3PN order, the PN approximation seems to 
present a convergent-like behavior (the region of validity 
increases with PN order). Similarly, PN approximants of 
order higher than 3PN present a divergent-like behavior 
(the region of validity decreases with PN order), proba- 
bly associated with the appearance of logarithmic terms 
in the series. 

Another result of this paper is the minimum number 
of multipoles that need to be retained to perform any 
type of comparison with PN theory (Table [J). We find 
that the number of multipoles depends on the velocity 
regime where the comparison takes place. In particu- 
lar, we find that more multipoles are required for low 
velocities. At 3.5PN order, the inclusion of the first five 
multipoles suffices to cover the entire velocity range. It 
will be interesting to explore the number of multipolar 
components that must be included in a numerical simu- 
lation of comparable mass binaries to get some required 
accuracy in the phasing of the waves. 

Although the method is generic and gauge- 
independent, the conclusions derived from the analysis 
of EMR systems are not necessarily generalizable to 
comparable mass systems. The number of multipoles 
needed to match the flux in the comparable mass case 
is probably smaller than suggested here. Moreover, the 
PN series is suspected to be much less accurate for EMR 
systems, which suggests that perhaps for comparable 
mass systems the onset of the divergent-like behavior 
observed here could occur at higher PN orders. 

Various extensions of this work should be possible in 
the near future. In particular, we would like to consider: 

1) Other physical systems: Although we have here 
studied only non-spinning EMR systems, the pro- 
posed method is very general and it can be applied 



15 



to more complex physical scenarios. For example, 
by modeling a small particle in orbit around a spin- 
ning (Kerr) black hole, one could study whether the 
black hole's angular momentum substantial affects 
the region of validity of PN theory. Other interest- 
ing scenarios would be comparable-mass, eccentric 
and spinning binaries. These systems, however, are 
much more difficult to model, both numerically and 
in PN theory. In particular, an accurate estimate 
of numerical errors in present numerical relativistic 
simulations is quite challenging (3f| [H, H3] • 

2) Different observables: We have here considered 
only the energy flux, because it is a well-defined, 
time-dependent observable. We could, however, 
also apply this method to the phasing of gravita- 
tional waves or the linear and angular momentum 
fluxes. One expects the asymptotic properties of 
PN theory to be independent of the quantity an- 
alyzed to determine it. If this were the case, the 
regions of validity found by analyzing, for example, 
the momentum flux or the phasing should be com- 
parable to those found here. Although unlikely, if 
this were not the case, then different PN-cxpanded 
quantities of the same order would possess different 
regions of validity, casting doubts on the asymp- 
totic structure and consistency of the PN approxi- 
mation. 

3) Different PN flavors: The method considered here 
does not depend on the PN approximant we use. 
We would like to study the region of validity of 
other (non- Taylor) PN approximants, such as the 
effective-one-body approach [p, B, Pade 0,111 and 
Chebyshev resummations [§.Tg. Ilfj||. In the EMR 
case, a fully convergent expansion of the flux might 
in principle be obtained via the Mano-Takasugi 
method [H, [H, 0, Sl|. To our knowledge, how- 
ever, the adiabatic energy flux has not been worked 
out explicitly in this approach, and their method 
does not seem to be straightforwardly generalizable 
to comparable masses. 

The studies suggested above would answer a number 
of questions and shed light on the structure of the PN 
series. Point (1) would reveal how the region of validity 
depends on the initial properties of the system. Based on 
previous studies [H, [l9| , one would expect this region to 
increase as the mass ratio approaches unity. Moreover, 
point (1) would shed light on the possible overlap of dif- 
ferent approximation schemes, such as PN theory, black 
hole perturbation theory and the close-limit approxima- 
tion [42j , providing support to asymptotically matched 
global solutions [43|, Sj, |45j. Point (3) would allow us 
to make relative statements about the regions of validity 
of different PN flavors. One could then verify whether 
different resummation techniques truly increase the re- 
gion of validity of the PN approximation (see e. g. [HI 
for criticism of this idea). In any case, a detailed study 



of the analytic structure of the PN approximation should 
greatly benefit our understanding of the dynamics of in- 
spiraling compact binaries. 
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APPENDIX: LOGARITHMIC TERMS 

The presence of logarithmic terms in the PN expansion 
[see e. g. Eq. (|17p] seems inherent to the PN approxima- 
tion in harmonic coordinates. The question we wish to 
answer is whether these terms affect the use of standard 
asymptotic techniques. The answer to this question is in 
the negative, and we explain why below. 

First, Arun et al. have argued that the logarithmic 
terms can be removed by a coordinate transformation for 
eccentric comparable-mass binaries. The new, so-called 
modified harmonic coordinates differ from standard har- 
monic coordinates to 1PN order and are regular every- 
where. Therefore, since the convergence or divergence 
properties of the approximation should not depend on the 
coordinate system used, it should not matter whether we 
do the analysis in harmonic or modified harmonic coor- 
dinates. In other words, if we were to repeat the analysis 
in modified harmonic coordinates, we should arrive at 
similar conclusions as those found here. 

Second, we can re-write the logarithmic terms in the 
PN approximation in a manner which clearly exhibits its 
asymptotic structure. Let us define v = lav, such that 
v € [-co, 0] as v € [0, 1]. Then we can write 

oo oo 

J2a P vHv p ) = e u ^b p v, (22) 

p— 6 p— 6 

where we have absorbed a factor of e p into b p . We then 
see that in the new variables this series is reminiscent 
of a standard asymptotic expansion with an exponential 
controlling factor (see e. g. |20l|). 

Third, we can straightforwardly show that the series 
is indeed asymptotic in some small velocity region. Let 
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us first note that any series w = Yin a « w ™ i s indeed 
an asymptotic series for small v, which is simply a re- 
statement that any Taylor expansion is an asymptotic 
series 20]. Similarly, it follows that a series of the 
form wi = Y n bn ln(w) is not an asymptotic series, be- 
cause ln('u) diverges as v — > 0. Fortunately, the PN 
expansion does not contain isolated logarithmic contri- 
butions. All logarithmic terms are multiplied by some 
power of the velocity, i. e. the expansion has the form 
w 2 = J2 n bnV n \n(v). As v — > 0, such a series indeed 
tends to zero and is well-behaved. 

Not only should the series be well-behaved as v — > 0, 
but it should also decay at the right rate, i. e. in such a 
way so that Eq. §5§ holds. Let us then look at the N = 6 
term in the PN expansion of the energy flux. Since we 



know that ciqv 6 decays at the right rate, we can compare 
b^v 6 ln(i>) to this quantity. One can easily verify that 
q,6V 6 is larger than the logarithmic term for v > 0.033, 
while the opposite is true for smaller velocities. When 
v ~ 10~ 4 and smaller, the logarithmic term dominates 
over the polynomial term, but such small velocities are 
excluded from our analysis anyway, since we only com- 
pute the numerical flux down to v = 10 -2 . Therefore the 
N = 6 term in the PN expansion of the energy flux is 
not only well-behaved as v — > 0, but it indeed decays at 
the right rate in the velocity region considered here. 

Based on these arguments, we conclude that there ex- 
ists a small- velocity near zone region, where logarithmic 
terms in the PN series do not affect the asymptotic anal- 
ysis described in this paper. 
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